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Abstract 

This paper demonstrates continuous optimization on the 
differentiable manifold formed by continuous constraint functions. 
The first order tensor geodesic differential equation is solved on the 
manifold in both numerical and closed analytic form for simple 
nonlinear programs. Advantages and disadvantages with respect to 
conventional optimization techniques are discussed. 

Introduction 


For many years operations research practitioners have sporadically 
considered using differential equations for optimization. Tanabe 
[1,2] took this approach into the realm of modern mathematics by 
deriving the properties of the tensor first order differential 
equation which can be used for optimization. This paper uses his 
theoretical foundations, extends his equation to a geodesic form, and 
demonstrates that nonlinear optimization problems can be solved 
both numerically and analytically with differential equations on the 
continuous differentiable manifolds formed by constraining a vector 
of continuous functions. 

To be consistent with the illuminating tensor notation of Gerretsen 
[3], row vector notation is used throughout this paper. Thus, a vector 
v with components a j in basis {b ^ , ... ,b n } may be expressed as 
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v = a B 


where 

a= [ai-..a n ] and B = 


All computer solutions were performed on a Macintosh Plus ™ using 
Lightspeed Pascal ™. 

Summary of Theory 

The following theory is based on concepts from differentiable 
manifolds and differential geometry as described by Thorpe [4] and 
Boothby [5], 

An n-dimensional manifold is a connected, locally compact space 
with a countable basis, each point of which has a neighborhood 
homeomorphic to euclidian n-space. 

A C k differentiable manifold is a manifold with additional 
mathematical properties imposed which permit the definition of 
compatible coordinate systems on the manifold which are mapped by 
diffeomorphic functions from the manifold into euclidian n-space. 

This mapping is depicted by Figure 1. 
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Figure 1 

Manifold Patch with Mapping into R 
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M c = { P I g(p) = c } can be shown to be a C 


k 


differentiable manifold. 


Such a manifold may be called a "constraint manifold" since it is 
totally defined given a set of C k constraint functions. 

The nonlinear programming problem may be stated as 


extremalize f(x) 
over x 

subject to g(x) = constant. 

For C k functions f and g with k > 0 and for the Jacobian matrix 


5 X 9 


~d^Q^...d n Q^ 
-5l 9nv--5n9m 


3g j 

where 3^ = ^ , 

there exists tensors N p and T p at the point p defined by 

Np = 5 T p g (5 p g 5 T p g)- 1 5 p g and T p =I-N p 

where 3 T p g is the transpose of 3 p g and the subscript p represents 
evaluation of the functions at p. 

Operations research analysts will recognize T p as a generalization 
of the projection matrix in the Rosen gradient projection algorithm 
[6] to evaluation at the tensor operating point p. 

Geometrically, T p projects onto the tangent space to the constraint 
manifold at p, and N p projects onto the normal space to the 
constraint manifold at p. 
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at the Tensor Operating Point p. 

A manifold M c is covered by a vector field if and only if at each point 
p there is a vector v(p). If f(x) is C k function over M c then the 

gradient 3 p f forms a C k vector field over M c . If f(x) and g(x) are C k 
functions over M c then 3 p f T p forms a C k vector field over M c . 

It is important to note that 3 p f T p is the restriction of 3 p f to M c . 
That means that 3 p fT p is in the tangent space to M c with origin at p. 

Tanabe [2] has shown that under appropriate second order conditions, 
the flow 

P = ± V T P 

extremalizes f(x) on M c by converging to a local extremum at steady 
state. 

A flow is a geodesic on a manifold M c if and only if all acceleration 
is normal to the manifold. 
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Consider the unit velocity flow 
±d P fT p 

P = ; • 

v3pfTp3 T p f 

on M c , then 

(a) p p T = 1, 

(b) p = 0, 

(c) any acceleration must be normal to M c , 

(d) the unit velocity flow is a geodesic on M c , 

(e) the unit velocity flow is the shortest distance on M c 

between any initial condition on M c and an extremum of 
M c , and 

(f) the distance between the initial condition and an 
extremum is the arc length of the path and hence a 
straight line on M c . 

A trivial extension of Tanabe's results [2] shows that under 
appropriate second order conditions the unit velocity flow 

±5pfT p 

P = , — — — 

Va p fT p a T pf 

extremalizes f(x) on M c by converging to a local extremum at steady 
state. 

At the extrema 

(a) p = 0, 

(b) 3 p f = 3 p fN p => All objective imparted velocity is 
normal to M c . 
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(c) If no two 3 p gj are colinear at p 

(1) N p has m eigenvalues Xj(p) such that A,,(p) = 1, 

(2) 3 p f can be expressed in the m dimensional basis of 
the eigenvectors associated with the m unit 
eigenvalues, 

(3) N p has n-m eigenvalues A,j(p) such that A,j(p) = 0, and 

(d) Nonlinear equations exist, which, if solvable, lead to a 
closed form solution. 

The Lagrange multipliers p, are 

m-(p) = V ^ T P g (^ P g 5 T P g) _1 

which gives 


P = d p f - P(P) ^ P g- 

At any point of the trajectory, the Lagrange multipliers (LLj are the 

normal components of the velocity which have been subtracted from 
the velocity imparted by the objective f(x) in order to restrict the 
trajectory to M c . 

Circular Manifold Analytical Solution 

Simple problems can be solved analytically for the dynamic path or 
just for the optimum points with this approach. The circular 
constraint manifold is solved below for a linear objective function 
to obtain the optimum point *p and optimum objective value f(*p). The 


problem is stated as 


maximize 

f(x) = a -) x-i + a 2 

over x-| and x 2 


subject to 

g(x) = x ! 2 + x 2 2 


6 



Evaluated at the point p we have 


3 p g = [ 2 Pl 2 p]> 

V = M 


N p = 3' p 9 ( a p 9 ^pg )' 1 a p 9 =2 


Pi P1P2 
P2P1P2 2 - 


Tp - I ' N p - r 2 


P2 -P1P2 
-P2P1 Pi 2 - 


and 


T p =^ 2 [ a 1 p 2 2 -a 2 Pi P2 _a i P2P1 +a 2$v‘ 


Setting 


V t p = ° 


we have 


a i P2 " a 2 Pi - 0 . 

The constraint provides 

* 2*2 2 

Pi + P2 = re- 
combining the two equations we have 


Pi = 


a 2 r 




2 +a 2 2 


a i r 

p 2 = and 

\ a-| 2 +a 2 2 

r/ * , 2 a-| a 2 r 

f(P) = , 

\ a-| 2 +a 2 2 

Circular Manifold Numerical Solution 
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The steepest descent trajectory on the circular constraint manifold 
is found by solving 

p = 3 p f T p 

= Jz [ a 1 P 2 2 ' a 2 Pl P 2' a 1 P 2 P 1 + a 2 f | ) 1 2 
The velocity along the manifold is 

V(p) = (P P T ) 1,z 

= (3pf Tp T p t 3pf T ) 1,2 
= (3 p f T p 3 p f T ) 1,2 . 

For an n dimensional sphere and hyperplane the velocity can be 
shown to be 


v(p) = ( a a T - |. 2 P i 1/2 ■ 

Using a fourth order Runga Kutta algorithm from Flanders [7] 
extended to multidimensional systems we have the following 
numerical trajectory. The radius has been set at 1, the components a 

set to 1, t is the trajectory measure, r 2 is the square of the radius 
of the sphere, v is the trajectory velocity at p(t), and p 1 , p 2 are the 
components of p(t). 

Note that the radius stays very close to 1.000000 throughout the 
trajectory. Consequently the algorithm follows the manifold 
closely. Also, since the unit velocity trajectory is not taken above, 
the velocity drops to zero as the optimum value is approached. This 
is actually desirable with the algorithm used. With a unit velocity 

trajectory it is difficult near the optimum to converge to the 
optimum because the constant distance step oversteps the optimum 
and hangs up. 


Table 1 

Steepest Descent of a Linear Objective 
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on a Circular Constraint Manifold of Radius One 


t 

r 2 

V 

P 

1 P 2 

0.00 

1.000000 

1.400000 

0.800000 

0.600000 

0.40 

0.999792 

1.296556 

0.365864 

0.930556 

0.80 

0.999000 

0.924766 

0.072559 

0.996863 

1.20 

0.998718 

0.572417 

0.360155 

0.932205 

1.60 

0.998855 

0.334783 

0.519319 

0.853910 

2.00 

0.999054 

0.191984 

0.604283 

0.7961 76 

2.40 

0.999294 

0.109368 

0.650076 

0.759405 

2.80 

0.999437 

0 . 0621 78 

0.675144 

0.737304 

3.20 

0.999584 

0.035324 

0.689081 

0.724397 

3.60 

0.999678 

0.020064 

0.696892 

0.716952 

4.00 

0.999755 

0.01 1395 

0.701300 

0.71 2694 

4.40 

0.999816 

0.006472 

0.703799 

0.71 0270 

4.80 

0.999851 

0.003676 

0.705214 

0.708889 

5.20 

0.999890 

0.002087 

0.706023 

0.7081 1 1 

5.60 

0.999912 

0.001 1 85 

0.706483 

0.707668 

6.00 

0.999930 

0.000673 

0.706745 

0.70741 8 

6.40 

0.999944 

0.000382 

0.706896 

0.707278 

6.80 

0.999952 

0.000217 

0.706981 

0.7071 98 

7.20 

0.999956 

0.0001 23 

0.707030 

0.7071 53 

7.60 

0.999959 

0.000070 

0.707057 

0.7071 27 

8.00 

0.999960 

0.000040 

0.707073 

0.7071 1 3 

8.40 

0.999961 

0.000022 

0.707082 

0.7071 04 

8.80 

0.999962 

0.000013 

0.707087 

0.7071 00 

9.20 

0.999962 

0.000007 

0.707090 

0.707097 

9.60 

0.999962 

0.000004 

0.707091 

0.707095 

10.00 

0.999962 

0.000002 

0.707092 

0.707095 

The unit velocity geodesic trajectory 

on the 

circular constraint 


manifold is found by solving 

b ml. 

p - v(p) • 

Spherical Manifold Numerical Solution 

The projection equations for the spherical constraint manifold are 
derived below. 
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The problem is stated as 

maximize f(x) = a ^ x-| + a 2 x 2 + a 3 x 3 

over x-| , x 2 and x 3 

subject to g(x) = x t 2 + x 2 2 + x 3 2 = r 2 . 

Evaluated at the point p we have 


3 p g = [2p 1 2p 2 2p] 3 


3pf = ^1 a 2 a i] 


N 


p " 


^ T P g (^ P g ^pg )' 1 ^ P g 



p 1 2 p 1 p 2 p 1 p 3 
P2PiP2 2 P2Ps 
-P3P1 P3P1P3 2 - 


T 


P 


I - N 



P 2 2 +P 3 2 -PlP 2 -PlP 3 
-P 2 Pl Pl 2 + P 3 2 -P 2 P 3 
- -P3P1 -P3P2 Pi 2 +P2- 


and 


d p f T p -^ 2 [a-)(p 2 2 +p 3 2 )-a 2 p-| p 2 -a 3 p-i p 3 , 

-a 1 P 2 P 1 +& 2 (p 1 2 +P3 2 )' a 3P2P3> 

' a 1 P3P1 ' a 2 P 3 P 2 +a 3 (P 1 2 + P2 2 )] 

The steepest descent trajectory on the spherical constraint 
manifold is found by solving 

P= V T P 

Using a fourth order Runga Kutta algorithm from Flanders [7] 
extended to multidimensional systems we have the following 
numerical trajectory. The radius has been set at 1, the components a i 

set to 1, t is the trajectory measure, r 2 is the square of the radius 
of the sphere, v is the trajectory velocity at p(t), and p 1 , p 2 , and p 3 
are the components of p(t). 


10 



Table 2 

Steepest Descent of a Linear Objective 
on a Spherical Constraint Manifold of Radius One 


t 

r 2 

V 

Pi 

P2 

P3 

0.00 

1.000000 

1.720465 

- 0.800000 

0.600000 

0.000000 

0.40 

0.999827 

1 .479468 

- 0.330427 

0.873341 

0.357660 

0.80 

0.998985 

0.902703 

0.107880 

0.842055 

0.527532 

1.20 

0.998620 

0.477765 

0.351083 

0.739594 

0.573028 

1.60 

0.998678 

0.242504 

0.468021 

0.665234 

0.580602 

2.00 

0.998835 

0.121738 

0.523763 

0.622780 

0.580217 

2.40 

0.998945 

0.060952 

0.550750 

0.600333 

0.578982 

2.80 

0.999054 

0.030496 

0.564019 

0.588830 

0.5781 14 

3.20 

0.9991 16 

0.015256 

0.570590 

0.583004 

0.577624 

3.60 

0.9991 73 

0.007631 

0.573867 

0.580078 

0.577372 

4.00 

0.999201 

0.00381 7 

0.575500 

0.578608 

0.577246 

4.40 

0.999230 

0.001909 

0.576321 

0.577875 

0.5771 87 

4.80 

0.999249 

0.000954 

0.576732 

0.577509 

0.577159 

5.20 

0.999263 

0.000477 

0.576938 

0.577326 

0.577148 

5.60 

0.999270 

0.000238 

0.577040 

0.577234 

0.577144 

6.00 

0.999273 

0.0001 1 9 

0.577091 

0.577188 

0.577142 

6.40 

0.999275 

0.000059 

0.577116 

0.577165 

0.577141 

6.80 

0.999275 

0.000030 

0.577129 

0.577153 

0.577141 

7.20 

0.999276 

0 . 00001 5 

0.577135 

0.577147 

0.577141 

7.60 

0.999276 

0.000007 

0.577138 

0.577144 

0.577141 

8.00 

0.999276 

0.000004 

0.577140 

0.577143 

0.577141 

8.40 

0.999276 

0.000002 

0.577141 

0.577142 

0.577141 

8.80 

0.999276 

0.000001 

0.577141 

0.577142 

0.577141 

9.20 

0.999276 

0.000000 

0.577141 

0.577141 

0.577141 


Since the unit velocity trajectory is not taken above, the velocity 
drops toward zero as the optimum value is approached. 


The unit velocity geodesic trajectory on the spherical constraint 
manifold is found by solving 

b Ip 

P ' v ( p ) ■ 

Table 3 illustrates the same solution as in Table 2 except the 
trajectory is along the unit velocity geodesic. A constant trajectory 
parameter step fourth order multidimensional Runga Kutta algorithm 
from Chua and Lin [8] was used. The trajectory parameter here was 
set at .01 and only every tenth step is shown until near the optimum 
value. After t=1.6 every step is shown. The optimum is between 
t=1.68 and t=1 .69, but the algorithm failed because there is no 



variable step end game to permit convergence to the optimum. Note 
the improved tracking of the manifold surface with r 2 > 0.999999. 
Note also that the trajectory parameter, in this case the arc length 
or minimum distance trajectory, reaches the optimum at a much 
smaller value corresponding to the shorter distance along the 
trajectory. 


Table 3 

Geodesic Descent of a Linear Objective 
on a Spherical Constraint Manifold of Radius One 


t 

r 2 

P 1 

P 

2 

0.00 

1.000000 

- 0.800000 

0.600000 

0.000000 

0.10 

1.000000 

- 0.747261 

0.661993 

0.058027 

0.20 

1.000000 

- 0.687055 

0.717371 

0.115474 

0.30 

1.000000 

- 0.619984 

0.765582 

0.171768 

0.40 

1.000000 

- 0.546719 

0.806143 

0.226345 

0.50 

1.000000 

- 0.467991 

0.838649 

0.278660 

0.60 

1.000000 

- 0.384587 

0.862776 

0.328192 

0.70 

1.000000 

- 0.297341 

0.878282 

0.374444 

0.80 

1.000000 

- 0.207123 

0.885013 

0.416955 

0.90 

1.000000 

- 0.114836 

0.882901 

0.455299 

1.00 

1.000000 

- 0.021402 

0.871968 

0.489095 

1.20 

0.999999 

0.165172 

0.824160 

0.541736 

1.30 

0.999999 

0.256448 

0.787763 

0.560056 

1.40 

0.999999 

0.345161 

0.743495 

0.572781 

1.50 

0.999999 

0.430426 

0.691799 

0.579782 

1.60 

0.999999 

0.511388 

0.633192 

0.580990 

1.61 

0.999999 

0.519215 

0.626973 

0.580791 

1.62 

0.999999 

0.526991 

0.620692 

0.580535 

1.63 

0.999999 

0.534714 

0.614349 

0.580220 

1.64 

0.999999 

0.542383 

0.607945 

0.579847 

1.65 

0.999999 

0.549996 

0.601481 

0.579416 

1.66 

0.999999 

0.557552 

0.594960 

0.578928 

1.67 

0.999999 

0.565043 

0.588387 

0.578382 

1.68 

0.999999 

0.572355 

0.581866 

0.577790 


Conclusions 

From the viewpoint of continuous optimization on constraint 
manifolds the theory and application of the nonlinear program 
reduces to the study of the zeros and the geodesic flows of the 
vector field 

dpfT p 

VdpfTpdV’ 
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The theory and examples in this paper provide the jumping off point 
for many additional studies. Both numerical and analytical dynamic 
solutions should provide valuable information about the geometry of 
the constraint manifolds by studying the natural straight lines on 
the manifold and the sensitivity of these geodesics to initial 
conditions. Many other possibilities exist. 

Preliminary numerical solutions indicate that without special 
tailoring of the differential equation solution algorithms, 
conventional minimum seeking optimization algorithms are faster. 
However, special algorithmic end games may reduce the advantage of 
conventional algorithms substantially. 

On the other hand, the amount of geometrical understanding of the 
constraint manifolds provided by the geodesic flows is far superior. 

Finally, as opposed to conventional nonlinear programming 
algorithmic approaches, the geodesic flow approach is in a 
mathematical form readily analyzed by recent advances in global 
dynamic theory. 
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